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Abstract 

This  report  is  a  description  of  the  NAVSPASUR  Operations  Center  process  of  re¬ 
ducing  the  receiver  station  data  to  local  direction  cosines  for  the  (new)  St.  Andrew’s 
cross  stations,  a  configuration  for  which  most  of  the  antennas  have  their  phase  cen¬ 
ters  aligned  along  north-east  and  south-east  directions.  The  following  observational 
parameters  are  also  determined:  d.  ection  cosine  rates,  Doppler,  and  chirp 

Most  of  the  main  steps  involved  are  the  same  as  for  the  older  cruciform  stations, 
whose  antenna  phase  centers  are  aligned  along  north-south  and  east-west  axes.  This 
algorithm  has  recently  been  modified  by  the  Analysis  and  Software  Department  of 
NAVSPASUR  for  use  with  the  St.  Andrew’s  cross  configuration.  The  direction  cosine 
algorithm  uses  a  walkup  and  least  squares  as  before,  but  differs  due  to  the  antenna 
configuration.  Those  modifications  were  in  effect  by  January  1,  1990,  the  freeze  date 
for  this  report.  It  is  hoped  that  the  present  report  will  elucidate  the  NAVSPASUR 
direction  cosine  processing  for  the  surveillance  community.  A  detailed  error  analysis 
of  the  steps  involved  in  the  direction  cosine  determination  processing  is  in  progress 
and  will  be  issued  eventually.  Most  of  the  present  report  will  be  incorporated  into  a 
longer  one  describing  the  NAVSPASUR  Operations  Center  process  of  satellite  element 
updates. 
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Introduction 


The  NAVSPASUR  sensor  is  a  multistatic  continuous  wave  radar  system  with  nine  stations 
(three  transmitters  and  six  receivers)  located  alohg  a  great  circle  path,  inclined  33.57  degrees 
to  the  equator  across  the  southern  United  States#(see  Fig.  1).  The  main  transmitter  is  at  the 
Kickapoo  Lake  Station,  and  the  wing  transmitters  are  at  Wetumpka  and  Gila  River.  The 
transmitters  broadcast  fan-shaped  beams  at  216.98  MHz,  and  the  unsynchronized  beams 
combine  to  form  the  NAVSPASUR  fence ,  which  is  very  narrow  in  the  north-south  direction, 
but  very  wide  in  the  east-west  direction.  When  an  object  such  as  a  satellite  enters  the 
fence  a  small  fraction  of  the  transmitted  radio  energy  is  reflected  to  one  or  more  of  the 
receiver  sites.  Interferometric  techniques  are  applied  to  determine  the  angles  of  arrival  of 
the  reflected  signal  at  each  receiving  station. 


Figure  1:  Map  of  NAVSPASUR  facilities 

The  data  from  all  receivers  are  collected  in  real  time  and  fed  via  dedicated  telephone  lines 
directly  into  the  central  data  processing  system  at  the  NAVSPASUR  Operations  Center  in 
Dahlgren,  Virginia.  At  the  Operations  Center  data  from  each  receiving  station  is  reduced 
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to  local  direction  cosines  for  each  object,  along  with  cosine  rates,  Doppler  and  chirp,  statis¬ 
tical  measures,  and  time  stamps.  The  NAVSPASUR  Operations  Center  maintains  satellite 
element  catalogs,  which  are  updated  with  the  reduced  data.  A  broad  overview  of  the  Oper¬ 
ations  Center  processing  is  given  in  the  flow  chart  of  Fig.  2.  The  main  part  of  this  report  is 
a  detailed  description  of  the  box  labeled  “RESOLUTION.” 


Figure  2:  NAVSPASUR  Operations  Center  data  processing 

The  local  coordinate  system,  called  East-North-Up,  is  defined  by  its  unit  coordinate 
directions.  Up  (2)  is  in  the  direction  of  the  local  vertical  as  established  with  a  plumb  line 
at  a  point  at  the  station,  and  the  antennas  are  in  a  perpendicular  plane,  established  with 
a  laser  survey.  East  and  North  coincide  with  the  directions  formed  by  the  antenna  phase 
centers  in  the  cruciform  array  shown  in  Fig.  3  (not  to  scale).  Moreover,  East  (1  or  E)  lies 
in  the  plane  of  the  great  circle  and  is  tangent  to  it;  North  ( y  or  N)  is  perpendicular  to 
the  great  circle.  If  a  satellite  is  in  the  direction  x,  from  the  center  of  the  ^ceiving  station, 
then  the  East  and  North  direction  cosines  for  the  satellite  are  defined  by  E  •  xs  and  N  •  xs, 
respectively.  These  quantities  are  the  cosines  of  the  angles  formed  by  xs  and  the  coordinate 
axes,  and  are  the  fundamental  observation  parameters  used. 
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Figure  3:  Typical  low  altitude  Cruciform  station  configuration  (not  to  scale) 

The  physical  deployment  of  the  antennas  is  intended  to  allow  for  the  unambiguous  deter¬ 
mination  of  direction  cosines  from  the  station  to  the  satellite.  The  antennas  are  configured 
as  an  array  of  collinear  half-wave  dipoles  and  are  positioned  to  form  interferometric  base¬ 
lines  of  varying  lengths.  For  each  antenna  the  radio  frequency  (RF)  signals  from  all  of  the 
dipoles  are  combined  coherently  (since  the  antennas  are  in  a  plane)  and  directed  through 
an  analog  filter  into  a  preamplifier.  The  filters  are  designed  to  remove  specific  signals  from 
local  transmissions  and  do  not  exist  at  all  sites. 

Four  of  the  receiving  stations  are  intended  for  detecting  low  altitude  satellites,  and  the 
other  two  are  for  high  altitude  satellites.  NAVSPASUR  has  recently  undertaken  a  program 
to  convert  the  antenna  array  configuration  from  the  cruciform  to  the  St.  Andrew’s  cross. 
The  low  altitude  St.  Andrew’s  cross  stations  are  illustrated  in  Fig.  4.  They  are  composed  of 
twelve  antennas  400  feet  long.  The  high  altitude  stations  are  composed  of  ten  antennas  each 
2400  feet  long.  The  antenna  configurations  for  low  and  high  altitude  differ  only  in  the  lack 
of  antennas  2  and  3  at  the  high  altitude  stations.  Antennas  1  and  4  overlap  by  1200  feet  in 
the  high  altitude  stations.  We  shall  use  the  name  Na  to  mean  the  number  of  antennas  at 
the  receiving  station  under  consideration  (ten  or  twelve). 

The  fundamental  quantity  measured  by  the  NAVSPASUR  sensor  system  is  the  set  of 
signal  phases  received  at  each  antenna  at  a  given  receiving  station.  The  raw  phases  from 
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Figure  4:  Low  altitude  St.  Andrew’s  cross  station  configuration 

each  of  the  antennas  are  combined  to  produce  phase  differences,  each  taken  with  respect 
to  the  designated  reference  antenna.  In  an  interferometer  system  these  relative  phases  are 
caused  by  differences  in  the  arrival  times  of  the  signal  at  the  different  antennas.  The  basic 
problem  is:  given  the  antenna  layout  and  the  measured  phases,  determine  the  satellite 
direction. 

For  a  given  satellite  pass  up  to  55  frames  of  signal  phases  are  recorded  at  a  54.93  Hz 
sample  rate,  allowing  the  determination  of  satellite  cosine  rates,  which  are  time  derivatives 
of  the  direction  cosines.  The  Doppler  and  chirp  of  the  received  signal  are  also  determined. 
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Each  receiver  station  has  a  primary  and  secondary  receiving  system  assuring  that  in 
the  event  of  equipment  fai'ure  the  station  is  still  able  to  perform  in  its  primary  mode. 
The  primary  mode  is  the  full  Doppler  processing  mode,  which  makes  nearly  all  of  the 
system’s  target  detections.  If  only  one  system  is  operational,  that  system  is  designated  the 
primary  system.  The  secondary  system,  when  available,  performs  in  the  half  and  quarter 
Doppler  processing  modes,  which  were  designed  to  detect  the  fainter  signals  of  higher  altitude 
satellites.  The  pertinent  details  of  the  full  Doppler  processing  mode  will  be  discussed.  The 
half  and  quarter  Doppler  processing  modes  and  the  secondary  system  as  a  whole  will  not  be 
treated. 

The  diagonal  arms  of  the  St.  Andrew’s  cross  are  at  45°  angles  to  the  east-west  direction. 
These  directions  are  referred  to  as  Lx  and  L2  f°r  northeast  and  southeast  respectively.  The 
direction  cosines  determined  with  respect  to  L\  and  Z,2  are  called  and  f2.  For  historical 
reasons  east-west  and  north-south  direction  cosines  are  desirable.  The  transformation  is 
accomplished  as  follows: 


m\  —  i  + 12) 

(1) 

where  rnx  and  m2  are  the  east-west  and  north-south  cosines. 

Other  observational  parameters  are  determined,  including  the  direction  cosine  rates, 
Doppler  and  chirp.  This  suite  of  processing  is  collectively  referred  to  as  direction  cosine 
processing.  An  overview  of  the  direction  cosine  processing  is  shown  in  the  flowchart  of  Fig.  5 
and  may  be  summarized  as  follows: 


•  Choose  an  observation  time  at  which  to  set  the  estimates. 

•  Smooth  the  phases  and  determine  the  direction  cosine  rates  with  least  squares. 

•  Determine  the  direction  cosines  with  walk  up  and  least  squares. 

•  Determine  the  Doppler  and  chirp  with  an  FFT  and  least  squares. 
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Figure  5:  Direction  cosine  processing  overview 


Section  1  is  an  overview  discussion  of  the  principles  of  interferometry  including  an  es¬ 
timate  showing  that  the  receiver  interferometer  near  field  effects  on  antenna  phase  do  not 
cause  significant  errors  in  the  direction  cosine  determinations.  The  format  of  the  phase  data 
which  arrives  at  NAVSPASUR  is  described  in  Section  2.  Details  on  the  observation  time  are 
in  Section  3.  Smoothing  the  phase  data  and  determining  cosine  rates  are  described  in  Sec¬ 
tion  4.  The  detailed  description  of  the  direction  cosine  processing  is  contained  in  Section  5. 
The  Doppler  and  chirp  calculation  is  described  in  Section  6.  An  appendix  summarizes  the 
least  squares  formalism.  Items  within  square  brackets  refer  to  bibliography  entries  at  the 
end. 
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Principles  of  Interferometry 


A  comprehensive  description  of  interferometry  may  be  found  in  [Kraus]  and  [Tho].  The 
simplest  interferometer  consists  of  a  pair  of  antennas  characterized  by  the  vector  from  the 
reference  antenna  to  the  remote  one.  This  vector  b  is  called  the  baseline.  In  the  limit  where 
the  source  of  the  radiation  is  infinitely  distant  (known  as  the  parallel  rays  assumption),  a 
plane-wave  signal  is  received  at  each  antenna  with  a  relative  phase  difference  which  depends 
on  the  direction  to  the  source  of  the  radiation.  This  is  illustrated  in  Fig.  6. 

/  / 
/  / 

/  / 


p  /  x  / 


Figure  6:  Two-element  interferometer 


The  angle  0  to  a  radiating  or  reflecting  satellite  is  given  by 

a  ~  b  p 
cos  V  =  s  •  y  =  7 
6  b 


(2) 


where  s  is  the  unit  vector  in  the  direction  of  the  satellite,  b  =  ]b|,  and  p  is  the  path  length 
difference  from  the  source  to  the  two  antennas.  If  the  received  voltages  are  expressed  as 
c0  and  tq  =  voe,2wi>,  then  under  the  assumption  of  parallel  rays  the  phase  difference  2ir4> 
between  the  received  voltages  at  the  antennas  is  2irp/\ ,  where  A  is  the  wavelength  of  the 
received  radiation.  In  terms  of  measured  quantities  the  direction  cosine  is  given  by 


cos  0  = 


(3) 


The  error  in  the  cosine  due  to  phase  measurement  error  is 

dcos  6  =  -rd<f>, 
b 
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indicating  that  longer  baselines  give  cosine  estimates  with  smaller  errors.  On  the  ocher  hand 
this  phase  difference  2tt4>  is  directly  measured  only  up  to  a  multiple  of  2tt.  As  a  consequence, 
cos  0  is  ambiguous  by  multiples  of  A/6.  This  is  known  as  the  cosine  ambiguity.  If  two  sources 
were  situated  with  direction  cosines  differing  by  A/6,  they  would  produce  identical  phase 
difference  data  on  the  antennas.  This  means  it  is  impossible  to  determine  from  phase 
measurements  whether  a  given  target  has  a  direction  cosine  iori  +  kX/b  for  k  an  integer. 
The  interferometer  is  also  said  to  have  grating  lobes  spaced  at  intervals  of  A/6,  although  the 
phrase  is  usually  reserved  for  ambiguities  in  the  multi-antenna  case.  Direct  phase  difference 
measurement  can  be  accomplished  if  the  baseline  is  shorter  than  one  wavelength. 

As  will  be  seen,  a  solution  to  the  competing  needs  for  short  and  long  baselines  is  based 
on  beginning  with  a  short  baseline,  and  applying  the  walkup  through  longer  baselines.  This 
overcomes  the  ambiguity  and  yields  the  intrinsic  accuracy  of  long  baselines. 

The  parallax  effect  due  to  the  satellite  being  in  the  near  field  of  the  interferometer 
receivers  does  not  cause  trouble  with  the  direction  cosine  estimates.  Parallax  arises  as  the 
error  in  the  measured  phases  due  to  the  parallel  rays  assumption  and  may  be  estimated  with 
trigonometry.  (See  Fig.  7).  Let  i?,  /?o,  and  Rx  be  the  slant  ranges  to  the  satellite  from  the 
center  of  the  interferometer  and  the  nearer  and  farther  antennas  respectively.  Then 

Rl  =  R2  +  62/4  +  bR  cos  6 

R\  =  R2  +  62/4  —  bR  cos  9 


so  that 


The  square  roots  may  be  approximated  by  an  application  of  the  binomial  expansion: 


(i  +  x)1/2  ~  i  +  \x  ~  ^x2  +  7£x3 

which  after  simplifying  and  discarding  terms  of  order  4  in  b/  R  yields 

1  b3 

Ro  —  Ri  ~  b  cos  6  —  -  —  cos  9  sin2  9. 

8  Rl 

The  estimate  for  cos#  is  Eq.  (2): 

- Ro-Ri  n  1  b2  .  2 

cos  9  -  - ; - ~  cos  9  —  -  —  cus  9  sin  9. 

b  8  R2 

The  random  variable  cos  9  sin2  9  has  a  mean  of  zero  and  a  standard  deviation  of  1  /2y/2.  This 
means  the  error  is  unbiased  and  has  a  standard  deviation  of  16^RJ ,  which  for  the  case  of  a 

baseline  of  1200  \/2  feet  and  a  satellite  range  of  100  km,  is  less  than  1.2- 10-6,  more  than  two 
orders  of  magnitude  smaller  than  the  nominal  rms  cosine  error  of  0.00025.  So  the  parallel 
rays  assumption  is  valid  for  all  St.  Andrew’s  cross  baselines  and  all  satellites. 
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2  Data  Format 


At  the  receving  stations  the  Interferometer  Subsystem  (IS)  measures  the  complex  voltages 
received  from  an  object  which  is  illuminated  by  one  of  the  transmitters.  Subsequent  data 
processing  (to  be  discussed)  determines  the  amplitudes  and  phases  of  the  signals. 

The  IS  includes  three  digital  filters  (DFs),  which  are  narrow  band  filters  passing  only 
the  selected  frequency  component  by  means  of  a  discrete  Fourier  transform  (DFT).  Tuning 
a  DF  is  equivalent  to  selecting  the  frequency  to  be  passed.  For  each  target  a  particular  DF 
is  tuned  to  the  Doppler  frequency  of  the  selected  target.  The  tuning  consists  of  storing  the 
appropriate  coefficients  for  multiplication  by  the  sample  values.  Each  DF  receives  all  of  the 
data  samples  from  all  of  the  antennas.  The  multiplexed  data  is  received  at  a  1.2  MHz  rate,  or 
75  kHz  for  a  single  antenna’s  data.  In  a  particular  DF,  sample  blocks  are  created  consisting 
of  4096  consecutive  samples  from  each  antenna.  Collecting  4096  samples  from  each  antenna 
requires  54.61  msec.  The  data  from  each  antenna  is  overlapped  three-to-one,  which  means 
each  sample  becomes  a  member  of  three  blocks,  and  one  block  for  each  antenna  is  completed 
every  18.2  msec.  This  is  illustrated  in  Fig.  8.  (Since  4096  is  not  evenly  divisible  by  three, 
a  new  block  is  started  after  each  of  two  consecutive  strings  of  1365  samples  and  another 
block  after  1366  samples  with  the  sequence  repeated.)  Each  block  is  windowed  with  a  four- 
term  Blackman-Harris  window,1  and  a  single  output  discrete  Fourier  transform  (DFT)  is 
performed  on  the  windowed  block.  A  single  complex  output,  the  voltage,  is  produced  every 
18.2  msec  for  each  antenna.  The  overlap  implies  that  the  voltages  have  correlated  error.  A 
maximum  of  55  data  samples  is  obtained  for  each  target. 


1365 

1365 

1366 

1365 

1365 

1366 

1365 

frame  1 - * 

i 

1 

1 

frame  2 


frame  3 


Figure  8:  Schematic  of  data  from  one  antenna 

The  noise  bandwidth  for  each  output  bin  of  a  DFT  performed  on  4096  samples  obtained 
at  a  75  kHz  sampling  rate  with  the  4-term  Blackman-Harris  window  is  36.6  Hz  with  bin 
centers  at  intervals  of  18.3  Hz,  so  that  at  worst  case,  a  satellite  could  be  9.15  Hz  away  from 
bin  center.  The  received  frequency  varies  during  a  satellite  pass  with  a  typical  Doppler  rate 

1  This  windowing  process  reduces  the  side-lobe  levels  to  —94  dB 
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of  —20  Hz/sec.  This  is  compensated  by  tuning  the  DF  one  bin  below  the  bin  in  which  the 
target  was  detected. 

Each  set  of  Na  complex  voltages  is  converted  to  Na  phases  and  one  representative  am¬ 
plitude.  This  comprises  one  frame  of  data.  The  full  set  of  frames  for  a  satellite  pass  is 
accompanied  by  a  header  frame,  a  set  of  integers  representing  the  relative  antenna  strengths. 
The  production  of  all  of  these  numbers  will  now  be  described. 

The  reported  amplitudes  and  antenna  strengths  are  calculated  as  follows:  Let 

s,(n)  =  |x|  +  |t/|  for  i=  l,...,Na 
and  n  =  1, . . . ,  55, 

where  x  and  y  are  the  real  and  imaginary  parts  of  the  complex  voltage  at  antenna  i  for  the 
nth  frame.  The  amplitude  for  the  nth  frame,  S(n),  is  the  Olympic  average2  of  the  amplitudes 
at  all  the  antennas  (the  s,(n)),  expressed  in  dB  above  1  mW  as 

ia(n)  =  -  [20  log  S(n)  -  SBIAS] , 

where  [■]  is  the  greatest  integer  function,  and  SBIAS  is  a  system  constant  for  scaling.  These 
integers,  za(n),  are  the  amplitudes  (typically  they  are  between  —160  and  —118  dBm).  With 
5  =  YlS(n),  and  C{  —  £ns,(n),  the  antenna  strengths  are  given  by 

C 

ianti  =  20  log  -£■  +  32 

limited  to  the  range  0  to  63.  The  ratio  C^/S  is  the  amplitude  ratio  for  the  zth  antenna,  and 
the  additive  constant  (32)  is  introduced  to  bias  the  numbers  away  from  zero. 

After  these  numbers  have  been  produced,  a  determination  is  made  of  how  many  frames 
to  send  to  the  Operations  Center:  the  numbers  S(n)  are  tested  in  reverse  order  until  a  value 
of  S(n)  is  above  the  threshold.  Let 

Nj  =  max{n  :  S(n)  >  Threshold}. 

The  frames  from  1  to  Nf  are  transmitted. 

To  accomodate  data  transmission  the  phase  at  each  antenna  is  quantized  to  6  bits,  as 
follows:  For  the  first  Nj  frames,  the  receiver  software  converts  (x,y),  the  real  and  imaginary 
parts  of  the  complex  voltage,  to  a  phase  (an  integer  in  [0,63])  using 

4>  =  [8—  -f  .5]  for  0  <  y  <  x. 
x 

2The  Olympic  average  of  a  set  of  numbers  is  obtained  by  deleting  the  largest  and  smallest  numbers  and 
taking  the  arithmetic  mean  of  the  remainder. 
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The  range  given  covers  only  the  first  octant;  the  second  octant  is  deduced  from 

tan'1  z  =  ff/2  —  tan-1  1  jz  for  z  >  1, 


which  yields 


<t>  =  16  -  [8-  +  .5] 

y 


for  0  <  x  <  y. 


The  rest  of  the  complex  plane  is  handled  with  simple  trigonometry.  The  error  for  the  noise- 
free  case,  e(4>)  =  <p  —  640  ( <p  in  rotations),  is  shown  as  a  function  of  tan-1  y/x  =  2ir<t>  for  the 
first  quadrant  in  the  plot  of  Fig.  9.  The  error  function  is  identical  in  the  other  quadrants. 


T  T 

0.0  */4  */2 

arc  tangent 

Figure  9:  The  phase  error  in  the  receiver  software  in  counts,  for  0  <  tan'1  y/x  <  7r/2. 

The  error  is  sometimes  larger  than  one  count  (one  count  is  5.635°),  even  without  noise. 

The  data  sent  to  the  ADDAS  at  the  Operations  Center  begins  with  a  header  line  contain¬ 
ing  the  receiver  doppler  bin  number,  a  set  of  integer  antenna  strengths,  and  a  time  stamp. 
Each  data  frame  consists  of  an  integer  representing  the  amplitude  and  twelve  integers  repre¬ 
senting  the  antenna  phases.  A  typical  data  record  is  included  as  Table  1.  The  ADDAS  is  the 
interface  between  the  telephone  lines  and  the  Operations  Center  computer  which  processes 
the  phase  data.  At  the  Operations  Center  the  integer  phase  is  divided  by  64,  expressing  the 
phase  in  rotations  6  [0,1).  We  will  use  the  notation  <j>i{n)  to  mean  the  phase  at  antenna 
i  for  the  nth  frame.  At  the  Operations  Center  the  ianti  are  converted  to  amplitude  ratios, 
wta,,  which  are  used  as  antenna  weights,  ,  * 

wta,  =  io  1'ant,“32. 

The  -3.2  removes  the  previously  applied  bias. 
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Table  1:  Typical  satellite  pass. 


-934  34  34  31  28  32  31  30 

39  14  47  19  61  61  18  36 

39  35  3  43  16  17  39  61 

-138  57  25  63  37  39  61  16 

-136  13  44  18  54  58  15  34 

-135  28  61  35  8  12  34  51 

-135  45  14  50  23  27  50  3 

-135  58  26  2  36  41  0  17 

-134  7  37  14  51  53  11  31 

-134  17  49  25  62  0  23  41 

-134  27  59  34  7  10  32  50 

-134  33  1  43  15  16  40  59 

-135  39  7  47  18  22  46  0 

-135  45  13  51  24  28  49  5 

-136  47  15  54  28  30  53  8 

-137  48  17  56  29  31  53  9 

-138  47  16  54  28  30  53  8 

-139  45  15  52  22  29  51  5 

-141  42  12  49  21  25  48  2 

-142  36  4  46  18  18  43  62 

-145  31  63  40  15  14  33  53 

-148  24  56  32  2  4  27  48 

-151  15  47  22  59  63  19  42 

-155  2  51  6  44  53  9  33 

-157  29  33  53  36  50  0  25 

-157  56  61  55  12  38  52  14 

-153  47  5  50  52  26  60  58 

-151  20  57  33  1  3  35  50 

-151  8  40  13  49  51  17  28 

-151  46  20  60  30  33  52  12 

-151  22  60  35  11  15  35  58 

-151  1  33  12  55  49  16  25 

-153  32  2  56  27  28  44  2 

-154  14  47  17  4  61  26  33 

-153  49  17  53  48  22  49  12 

-158  5  4S  21  45  12  21  28 

-156  8  16  8  26  29  40  1 

-155  10  47  5  6  35  37  48 


28  34  33  32  30  890404  900110.901 
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3  The  Observation  Time 


A  specific  frame  is  selected  at  which  to  base  the  direction  cosine  estimates.  It  is  desirable 
that  the  frame  (called  NFTO)  be  at  a  time  when  the  amplitude  of  the  signal  is  high,  and  that 
the  satellite  north-south  cosine  be  near  zero.  The  algorithm  for  the  selection  of  the  frame 
NFTO  follows:  Let  the  integer  NF1  be  the  number  of  the  frame  which  is  the  first  occurrence 
of  the  highest  amplitude.  Integer  NF2  is  a  weighted  mid-point  for  the  satellite  pass  which  is 
defined  to  be  the  maximum  M  satisfying 

M  t  Nf 

£  Wtfn  <  9  £  Wtf" 
n=l  Z  n=l 

where  Nj  is  the  number  of  frames  in  the  pass.  The  desired  frame,  NFTO,  is  the  truncated 
average  of  NF1  and  NF2.  The  observation  time  is  the  time  of  frame  NFTO. 

The  weight  function  will  be  used  in  the  phase  smoothing,  and  is  obtained  in  the  following 
manner:  let  ia  =  ia(n)  be  the  amplitude  for  the  nth  frame,  then  let 

x  =  (  10"  if  *<*>-160; 

l  1  otherwise 

(a;  is  inversely  proportional  to  the  amplitude  down  to  —160  dBm);  then  the  weight  wtfn  for 
the  nth  frame  is 

wtfn  =  — ? - : - - 

£LoC,x*  + 0.0045* 

where  the  coefficients  are  given  by 

co  =  0.000135 

d  =  0.016142 

c2  =  0.313793 

c3  =  -1.174551 
c4  =  1.357612 

c5  =  -0.374264, 

and  0.0045  is  1  /64th  of  1  / -\/T2,  the  standard  deviation  of  the  error  due  to  6  bit  quantization 
of  the  unit  interval,  [0, 1].  The  form  of  the  weight  function  was  determined  by  Mr.  Richard 
Smith  (private  communication),  the  head  of  the  NAVSPASUR  Analysis  and  Software  De¬ 
partment,  since  retired.  He  used  a  polynomial  fit  of  the  variance  vs.  amplitude  for  the  data 
taken  at  the  beginning  of  Phase  4  of  the  Modernized  Receiver  System.  The  weight  function 
is  shown  in  the  plot  of  Fig.  10. 
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-160  -120 
Figure  10:  The  weighting  function  of  frame  amplitude 

4  Cosine  Rates  and  Smoothed  Phases 

The  center  of  the  antenna  array  is  selected  as  the  reference  point.  An  iterative,  weighted 
least  squares  differential  correction  procedure  is  performed  on  the  phase  data  to  produce 
smoothed  phases  on  all  the  antennas  at  frame  NFTO,  and  for  all  the  frames  at  the  reference 
point  (but  see  the  next  paragraph.  The  latter  is  called  the  phase  profile.  The  differential 
correction  also  provides  satellite  direction  cosine  rates.  All  of  this  is  described  in  this  Section. 
The  smoothed  phases  at  the  time  of  frame  NFTO  are  used  to  calculate  the  direction  cosines 
in  the  next  Section.  The  use  of  the  phase  difference  profile  to  determine  the  Doppler  and 
chirp  is  discussed  in  Section  6. 

The  cosine  rates  and  smoothed  phases  are  determined  in  two  or  three  least  squares 
iterations.  As  will  be  seen,  the  phase  profile  is  not  directly  solved  for;  rather  the  differential 
correction  determines  the  phase  difference  profile,  the  phase  profile  minus  the  phase  at  NFTO. 
Usually  three  steps  are  performed,  but  during  heavy  traffic,  provision  is  made  to  limit  the 
processing  to  two  iteration  steps.  Also,  as  will  be  described,  if  after  the  second  iteration, 
the  estimates  are  good  enough,  there  is  no  third  iteration.  The  following  description  is 
summarized  as  the  flow  chart  of  Fig.  11. 

In  the  first  step  only  the  phases  at  NFTO,  <pi  =f  <^,(NFT0),  and  the  cosine  rates  are  adjusted; 
0.0  is  used  for  the  initial  values  of  the  cosine  rates.  The  measured  phases  at  all  Na  antennas 
for  all  Nf  frames  are  the  observations;  the  calculated  values  from  the  model  are  given  by 

c,(n)  =  fa  +  A<j>(n)  +  rhtXitn  +  m2y,tn , 

where 
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<t>i  =  0i(NFTO)  is  the  model  phase  at  antenna  i  at  the  observation  time, 
rhi  and  m2  are  the  east-west  and  north-south  cosine  rates, 

A 4>(n)  =  <f>re{ (n)  —  <pTe{  is  the  model  phase  difference  profile, 

4>rei{n)  is  the  phase  profile,  the  phase  at  the  reference  point  for  the  nth  frame, 
x,  and  yi  are  the  antenna  coordinates  with  respect  to  the  reference  point, 
and  tn  is  the  time  difference  between  the  nth  frame  and  NFTO. 

Only  the  first  two  are  adjusted  in  the  first  iteration. 


Figure  11:  Flow  chart  summary  of  phase  smoothing  and  cosine  rates. 
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The  term  least  squares  means  a  process  of  minimizing  a  squared  difference.  The  most 
common  example  is  fitting  a  curve  to  measurements,  but  there  are  many  scientific  and 
engineering  applications  of  the  technique,  which  may  be  viewed  as  a  generalization  of  the 
curve-fitting  problem.  The  theory  is  developed  in  many  places  including  pp.  123-5  of  [Bate] 
and  is  summarized  in  an  appendix. 


The  least  squares  correction  process  requires  the  solution  to  the  equation 


AiWA  A P  =  A'W(0  -  C), 


where  A  is  the  matrix  of  partial  derivatives  of  calculated  values  with  respect  to  the  param¬ 
eters: 

dc,(n)  dci(n) 


=  St 


d<f>:  11  3(mi,m2) 


=  (*;*», &*„)• 


A  is  given  by 


A  = 


/  Ik  xti  yt\  \ 

Ik  zb  ift  2 
^  Ik  ylNf  ) 


where  I^a  represents  the  identity  matrix  of  order  Na,  and  x  and  y  are  the  column  vectors 
of  the  x  and  y  coordinates  of  the  antennas.  The  weight  matrix  W  is  diagonal  with  entries 
equal  to  the  relative  weight  assigned  to  each  observation, 


W  =  diag  {wtu,...  ,wt^,i,wtlf2,...  ,wt^,2,.  ..,wt^x}  • 


Here  wtnj  is  the  weight  assigned  to  the  phase  measurement  at  antenna  i  for  frame  n,  and  is 
given  by  the  product,  wtfnwta,,  of  the  weight  for  frame  n  (determined  in  Section  3)  and  the 
weight  for  the  zth  antenna  (Section  2).  A P  is  the  vector  of  improvements  to  the  parameters 
P,  and  0  —  C  is  the  vector  of  residuals,  <f>i(n)  —  c,(n),  which  are  normalized  to  [—0.5, 0.5] 
by  adding  or  subtracting  1,  which  avoids  the  usual  difficulties  in  unwrapping  the  measured 
phases. 


The  matrix  AiWA  is  known  as  the  inverse  covariance  matrix  because  when  the  weights 
are  the  reciprocals  of  the  variances  of  the  observations,  ( AlWA )  1  gives  the  covariance 
matrix  for  the  new  estimates.  In  particular,  its  ith  diagonal  entry  is  the  variance  (square  of 
the  standard  deviation)  of  the  adjusted  estimate  of  the  ith  parameter,  P,  +  AP,. 


A‘Wi4  may  be  expressed  as  a  block  matrix 


A*WA 


m  n  m12 
mi2l  m22 
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where  mn  is  diagonal,  Na  x  Na,  with  element  ii  given  by  £n  wtni.  The  Na  x  2  matrix 

mu  =  ^wtnif„(x,-  ».-)> 

n 


and 

is  2  x  2. 


m22  =  EE  wtni<* 

*  n 


X?  X.y,  \ 

x.y,  V?  ) 


The  vector  AtVV(0  —  C )  has  +  2  elements 

£n  wt ni{<f>i(n)  -  c,(n)) 

£.  £n  Wtn.in  (<Mn)  ~  C,(n))  ^  ^ 

with  the  residuals,  (<?^(n)  —  c<(n)),  normalized  to  [—.5,  .5]. 

The  block  structure  of  AlWA  makes  the  matrix  inversion  fairly  straightforward,  resulting 
in  adjustments  to  be  made  to  <f>i  and  the  cosine  rates. 

The  second  iteration  begins  with  the  adjusted  values  as  well  as  the  unadjusted  history 
of  phase  differences  A(f>(n),  which  are  now  also  to  be  adjusted.  The  calculated  values  take 
the  same  form: 

c.'n)  =  <f>,  +  A<f>{n)  +  mix<fn  +  m2yttn. 

The  matrices  involved  are  considerably  larger  than  on  the  first  iteration,  since  there  are  Nj 
more  parameters.  Also,  only  observations  meeting  the  tolerance,  i.e.,  that 

wt„,  |0j(r»)  -  c,(n)|2  <  T0i, 

are  included  in  the  sums,  where 

r"l  =  4mlX{IVaJV/-4i-^-2'1} 

and  S2  =  £,  £„  wt„,  ( <t>i{n )  —  c,(n))2  is  the  weighted  sum  of  the  squared  residuals  from  the 
first  iteration.  The  matrix  A  now  includes  the  partial  derivatives 

dcj{n)  _ 
dA “  nm’ 

so  A  is  given  by 

lNa  Ex  xtx  yt , 

4va  E2  xt2  yt2 

iNa  En;  xtsf  yiN, 
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where  Ek  is  the  Na  x  Nj  matrix  whose  kth  column  is  all  Is,  and  the  other  entries  are  0. 

The  submatrix  mn  is  no  longer  diagonal,  and  now  of  size  Na  +  Nj  x  Na  +  Nj,  in  the  form 


mn  = 


Enwtr 


wtr 


(wtn.)'  E.  wtm  )  ' 

where  wtni  (unsummed)  represents  the  Na  x  Nj  block  of  individual  weights.  Also 

(  En  Wt nitn{x,  Vi)  \ 
mu  =  I  , 

\  Et  Wtm'in(;T,'  Vi)  J 

is  of  size  Na  +  Nj  x  2,  and  m22  is  as  before. 

The  vector  AlW{0  —  C )  has  Na  +  Nj  +  2  elements 

/  En  wtn,  (^.(n)  -  Ci(n)) 


E.  wt„,  {<f>i(n)  -  ct(n)) 
E,  En  wt  n,tn(<t>i(n)  ~  ct(n))  ( 


Vi  J 


At  the  Operations  Center  an  approximation  to  AlWA  is  made;  mu  is  taken  to  be  diagonal: 


mu  = 


En  wt„ 
Oni 


Om 

E<wt„ 


where  represents  the  k  x  j  matrix  of  zeroes. 

If  En  wtni  =  0  f°r  some  i,  meaning  none  of  the  phase  measured  at  antenn?  met  the 
tolerance,  or  E<  wtn«  =  0  for  some  n,  meaning  none  of  the  phases  measured  at  frame  n  met 
the  tolerance,  the  matrix  A'WA  as  given  would  be  singular.  Instead  the  system  of  equations 
is  purged  of  the  row  and  column  containing  the  term. 


A  third  iteration  is  made  if  traffic  has  not  been  too  heavy,  and  if  the  dot  product  of  the 
quantities  AtW(0  —  C)  and  the  applied  correction  is  more  than  20%  of  52,  the  weighted 
sum  of  the  squared  residuals  which  met  the  tolerance  during  the  second  iteration.  If  a  third 
iteration  is  to  be  made,  a  new  tolerance  is  computed: 

f  S2NfNa 
\  So  (So  —  Nj  —  \  —  2)  J 

where  So  is  the  number  of  observations  that  met  the  tolerance  in  the  second  iteration. 
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The  variances,  Var^,  of  the  estimates  of  d>,-  are  set  to  numbers  proportional  to  the  diagonal 
elements  of  (^‘WA)-1  and  recorded,  so  that  new  weights  for  the  antennas  can  be  used  for 
direction  cosine  estimation. 

The  resulting  phases,  fa,  are  calibrated  by  subtracting  the  calibration  phases.  The 
adjusted  (and  calibrated)  values  will  be  referred  to  with  the  same  symbols,  (f>,  for  the  Na 
phases  at  frame  NFTO  (whether  normalized  to  [0,1]  or  to  [—.5,  .5]),  m.\  and  m2  for  the 
direction  cosine  rates,  and  A  fan)  for  the  phase  difference  history  at  the  reference  antenna. 


5  Direction  Cosines 


The  walkup  technique  is  a  systematic  step-by-step  resolution  of  the  direction  cosine  from 
the  phase  differences,  beginning  with  an  estimate  due  to  the  phase  difference  on  the  shortest 
baseline  and  adjusting  the  estimate  in  steps  through  sequentially  longer  baselines.  The 
procedure  is  performed  on  both  arms  of  the  St.  Andrew’s  cross  to  determine  the  direction 
cosines  with  respect  to  each  of  the  arms.  In  this  Section  the  theory  is  explained  and  then  the 
implementation.  The  direction  cosines  thus  determined  are  improved  with  a  least  squares 
differential  correction  procedure. 


5.1  The  Walkup  Method 


This  algorithm  has  recently  been  modified  by  the  Analysis  and  Software  Department  [Bales] 
of  NAVSPASUR  for  use  with  the  St.  Andrew’s  cross  configuration.  Let  6,  =  6,/A  be  the  2th 
baseline  (in  wavelengths)  in  ascending  order  of  length,  and  let  <t>,  be  the  calibrated  phase 
difference  in  rotations  between  the  phases  measured  on  the  two  antennas  comprising  baseline 
6,.  In  the  Operations  Center  software,  the  phase  differences  are  normalized  to  [0, 1],  but  we 
shall  assume  <p ,  €  [—.5,  .5],  which  makes  the  description  easier.  The  initial  estimate  of  the 
direction  cosine  is  as  given  in  Eq.  (3):  Cj  =  d>i/6j.  The  iteration  step  of  the  walkup  is  as 
follows:  Assume  we  have  c,_t,  the  cosine  estimate  after  processing  baseline  i  —  1,  and  let 
$  =  c,_ j 6,  and  $  =  [$]  +  fa.  If  the  absolute  value  of  $  —  $  is  greater  than  0.5,  then  $ 
is  adjusted  toward  $  by  ±1.  The  quantity  $  is  the  total  phase  difference  on  baseline  6,, 
composed  of  an  integer,  [$],  and  fractional  part,  fa.  The  new  cosine  estimate  is  given  by 
c,  = 

In  order  that  cx  be  uniquely  determined  from  the  phase  difference  on  baseline  6,,  be.,  that 
there  be  no  cosine  ambiguity  at  the  first  step,  it  is  necessary  that  cf  =  (d>j±l)A/6j  represent 


20 


values  larger  than  1  / \/2  which  are  unphysical,  since  direction  cosines  on  the  diagonal  base¬ 
lines  are  between  approximately  45°  and  135°.  The  condition  is  met  provided  b\  <  X/\/2. 
Similarly,  in  order  to  avoid  an  ambiguity  at  each  step  in  the  walkup  procedure,  it  is  necessary 
that  the  sequence  of  baselines  (6<)  grow  slowly.  The  error  in  $  is  6,/6,_ t  (=  6, / 6^ _ x )  times 
the  error  in  the  phase  difference,  fa.  The  restriction  that  &,/&,_!  <  2  keeps  the  error  on  <J» 
low  enough,  it  is  hoped,  to  determine  the  integer  [$]  correctly  sufficiently  often. 


5.2  Implementation 

Physical  baselines  along  Lx  and  L2  do  not  satisfy  the  conditions:  the  shortest  is  20\/2  feet, 
and  there  are  ranges  where  we  would  have  6,  >  26;_!.  The  lists  of  baselines  for  Lx  and  L2 
are  in  the  first  columns  in  each  half  of  Table  2.  This  is  overcome  with  virtual  baselines, 
created  using  second  (and  more)  differences.  For  example,  there  are  baselines  along  L2  of 
20\/2  and  44\/2  feet,  whose  phase  differences  can  be  differenced  to  form  a  virtual  baseline 
of  24 \/2  feet. 

A  virtual  antenna,  VI,  is  introduced  to  fill  the  missing  corner,  at  the  point  (1200,  1200) 
(see  Fig.  4),  with  a  phase  (in  this  paragraph  we  are  using  if),  for  the  phase  at  antenna  i 
rather  than  the  phase  difference  on  baseline  b,.) 

0vi  =  fractional  part(<^i  +  fa2  —  <^4), 

and  variance  ^ 

VarVi  =  -  (Varx  +  Var12  +  Var4) , 

which  may  be  compared  with  the  theoretical  value 

Varvi  =  Vari  -f  Vari2  -f  Var4. 

The  Zj  baselines  ad>.  th  antenna  V\  are  shown  in  the  second  column  of  Table  2.  Another 
virtual  antenna,  V2,  is  introduced  at  the  center  of  the  X;  it  will  be  assigned  a  phase  equal 
to  the  average  of  the  phases  at  antennas  1  and  12,  i.e., 

<t>\i  —  ^  (<Ai  +  ^12)  or  <t>\i  —  2  (^1  +  ^12  +  1)  1  (4) 

and  variance 

Varv2  =  -  (Vari  +  Varx2) , 
which  may  be  compared  with  the  theoretical  value 

Varv2  =  7  (Vari  +  Var12) . 

4 
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Table  2:  St.  Andrew’s  cross  baselines 
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The  baselines  added  with  antenna  V2  are  shown  in  Table  2. 

As  indicated  4>v 2  is  ambiguous  by  ±0.5:  suppose  (f>x  =  —0.4  and  <t>  12  =  0.4;  then  one 
can  have  <f>v2  =  0.0,  but  <f>v2  —  ±0.5  is  equally  possible.  This  has  assumed  the  <fi,  are 
normalized  to  [—0.5, 0.5],  but  other  normalizations  are  analagous.  Actually  the  ix  walkup 
begins  without  a  value  assigned  to  <j>\2.  The  phase  differences  for  baselines  including  V2  are 
tentatively  chosen  based  on  whatever  value  <j>y 2  happens  to  have  (random).  Near  the  end  of 
the  ix  determination,  (j>\ 2  is  recomputed.  We  will  see  that  there  is  no  harm  in  this. 

A  total  of  42  baselines  are  used  in  the  walkup,  31  for  £x  and  12  for  i2.  Three  of  the  L2 
baselines  are  used  on  Lx;  these  are  20,  44,  and  64  feet  (times  \/2).  If  it  is  true  that  at  frame 
NFTO,  the  north-south  direction  cosine  is  small,  then  from  Eq.  (1)  we  must  have  that  £x  ~  t2. 
In  this  case,  the  phase  differences  which  would  be  produced  on  short  baselines  on  Lx  and 
L2  will  be  approximately  equal.  This  is  the  justification  for  using  short  L2  baselines  as  if 
they  were  along  L\.  As  will  be  seen  the  rather  longer  L\  baseline  of  194 \/2  is  used  for  the 
first  i2  walkup  step.  The  last  columns  in  each  half  of  Table  2  show  the  baselines  contributed 
through  the  use  of  Lx  and  L2  baselines  (and  virtual  baselines)  on  L2  and  Lx. 

The  algorithm  for  determining  the  direction  cosines  from  the  smoothed  phases  at  the 
reference  time  (frame  NFTO),  as  described  in  the  rest  of  this  Subsection,  is  summarized  as 
the  flow  chart  of  Fig.  12.  The  list  of  the  phase  differences  used  for  the  baselines  is  shown  in 
Table  3,  where  the  common  factor  of  >/2  has  been  suppressed.  The  shortest  baseline  used 
in  NAVSPASUR  direction  cosine  processing  is  2\/2  feet  which  is  sufficient  to  avoid  a  cosine 
ambiguity  at  the  first  step,  i.e.,  satisfies  the  condition  bx  <  X/y/2,  since  X  is  more  than  4 
feet.  With  this  set  of  baselines  on  Lx,  we  have  bi  <  26,_i  except  for  the  shortest  pair  of 
baselines. 


The  cosine  on  Lx,  £x,  is  determined  by  walkup  through  all  baselines  indicated  in  the  Table 
(The  ones  up  to  194 \/2  inclusive  are  independent  of  <j>v2 ),  testing  each  successive  estimate. 
If 

|c;  c»-i|  ^  -25/6,, 


then  a  weighted  least  squares  estimate  is  performed  (as  described  in  the  next  Subsection) 
on  the  baselines  and  phases  up  to  i  (but  not  including  the  ones  that  depend  on  <pvi)  to 
redetermine  c^,  using  c,_i  as  the  initial  estimate.  Since  </>v2  has  not  yet  been  chosen  (see 
above),  baselines  dependent  on  <f>v2,  c,  will  include  a  random  error,  which  is  likely  to  trigger 
the  least  squares  process  to  recalculate  c,  without  using  baselines  dependent  on  4>v i- 


After  the  walkup  has  proceeded  through  the  longest  L\  baseline,  regardless  of  whether 
the  resulting  c,  meets  the  0.25/6j  tolerance,  a  least  squares  fit  is  performed  for  ix  using 
baselines  independent  of  <f>v2-  This  allows  a  determination  of  <f>v2,  choosing  from  Eq.  (4) 
the  value  for  <^2  whose  phase  difference  with  respect  to  antenna  4  (see  Fig.  4)  most  closely 


23 


Table  3:  Antenna  pairs  used  to  create  baselines 


Lx 

Baseline 

Antennas 

l2 

Baseline 

Antennas 

1200 

Vl-4 

1200 

12-1 

1110 

9-4 

1180 

11-10 

1082 

8-4 

1136 

10-1 

916 

9-5 

912 

12-6 

888 

8-5 

892 

11-6 

830 

7-4 

848 

10-6 

636 

7-5 

600 

12-V2 

600 

V2-4 

580 

11-V2 

510 

9-V2 

536 

10- V2 

482 

8-V2 

312 

V2-6 

406 

V2-5 

288 

6-1 

370 

Vl-7 

194 

5-4 

280 

9-7 

64 

12-10 

252 

8-7 

44 

11-10 

230 

7-V2 

20 

12-11  | 

194 

5-4 

118 

Vl-8 

90 

Vl-9 

64 

12-10 

44 

11-10 

28 

9-8 

20 

12-11 

18 

11-10  Vl-9  12-10 

16 

11-10  9-8 

14 

12-10  Vl-9  12-11* 

12 

12-11  12-10  11-10  9-8 

10 

9-8  11-10  Vl-9  12-10 

8 

9-8  12-11 

6 

Vl-9  12-10  12-11 

4 

12-11  11-10  9-8 

2 

12-10  Vl-9  9-8 

*  The  20\/2  foot  baseline  is  used  twice. 


25 


matches  the  phase  difference  calculated  for  a  600\/2  foot  Lx  baseline.  The  phase  differences 
which  use  <f>v 2  are  recomputed.  Provision  is  made  for  a  final  least  squares  fit  to  be  performed 
on  ix,  this  time  using  all  the  baselines.  The  final  least  squares  fit  on  £x  is  not  done  if  the 
penultimate  estimate  is  considered  adequate.  See  the  next  Subsection. 

This  final  value  of  i\  is  used  as  the  starting  value  for  1 2 ,  on  which  a  walkup  is  performed, 
first  to  the  194 \/2  foot  baseline  of  Lx,  and  then  on  longer  L 2  baselines.  A  least  squares  fit 
is  performed  on  i2.  Finally,  these  values  are  rotated  back  into  north-south  and  east-west 
components  using  Eq.  (1). 


5.3  Least  Squares  Adjustment 


During  the  processing  described  in  the  previous  Subsection,  the  direction  cosines  ix  and  £2 
are  corrected  with  a  least  squares  fit  to  the  phase  differences  </>,•  using  fa  =  £,6,  (derived  from 
Eq.  (3))  as  the  model.  The  partial  derivatives  involved  are 


dfa  _  r 
at,  " 


Since  there  is  only  one  parameter  being  adjusted  (£3  and  £2  are  handled  separately)  the 
differential  corrections  to  £x  and  £2  are  the  solutions,  A£3  and  A£2,  of  scalar  equations  and 
given  by: 


A£j  = 


£wtb jbjjfa  -  bjlj) 
£  wtb,-6- 


where  j  =  1,2  and  the  sums  run  over  the  baselines  used. 


(5) 


There  is  no  tolerance  test  for  the  residuals  in  the  least  squares,  but  not  all  of  the  baselines 
are  used.  During  the  £x  walkup,  if  a  least  squares  fit  is  performed  for  Lx  baselines  up  to 
(and  including)  194  feet  (this  happens,  for  example  if  |q  — c,_j|  >  .25/6,),  all  baselines  up 
to  the  present  one  are  considered,  including  the  three  short  L2  baselines  and  all  the  virtual 
baselines  (from  2\/2  to  18\/2  feet).  For  walkup  on  longer  baselines,  the  shortest  nine  and 
the  L2  baselines  are  excluded,  so  it  starts  with  28\/2;  it  does  not  use  the  ones  that  use  <p\2. 
After  the  first  least  squares  for  the  1200\/2-foot  baseline  (which  does  not  include  the  L2  or 
(f>V2  baselines),  then  a  test  is  made.  If 


^  Ewtb j(fa  -  bt£j)2 

£wtbi6,2 


<  10 


-7 


(6) 


where  JV6  is  the  number  of  baselines  involved  in  the  least  squares  (14  in  this  case),  then 
the  new  value,  £x  +  A£j,  is  accepted,  and  there  is  no  final  least  squares  step  for  Cx.  If  the 
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inequality  does  not  hold,  then  after  the  final  least  squares  step  for  £lt  the  left  hand  side  of 
Eq.  (6)  is  recomputed.  The  numerator  is  the  sample  variance  of  the  previous  estimate  of  lXy 
and  the  denominator  is  the  theoretical  variance  of  the  new  estimate.  The  final  (extra)  least 
squares  correction  step  of  i\  (if  performed)  includes  the  baselines  which  use  <^V2»  it  having 
been  determined  in  the  meantime.  During  i 2 ,  including  its  final  least  squares  correction 
step,  if  the  walkup  step  requires  least  squares,  all  baselines  on  L 2  (up  to  the  present  one) 
are  used. 

The  weight,  wtb<,  for  the  baseline  6,  is  determined  from  the  variances,  Varal  and  Vara2, 
of  the  (real  or  virtual)  antennas,  al  and  a2,  comprising 


wtb,  =  1/  (Varal  -f-  Vara2) . 


The  virtual  baselines  are  assigned  weights  equal  to  the  sum  of  the  weights  of  the  baselines 
comprising  the  virtual  baselines,  except  that  the  baseline  of  length  12\/2  feet  gets  the  weight 

wt12  =  wt28  +  ^  (wt20  +  wt44  +  wt64 ) , 

where  the  subscripts  identify  the  baseline  length  (except  for  the  \/2).  The  14-\/2  feet  baseline 
includes  the  20\/2  feet  baseline  twice,  so  its  weight  is  added  twice  in  the  sum  for  the  weight 
of  the  virtual  baseline. 

For  the  weights  assigned  to  virtual  baselines  to  be  mathematically  consistent  with  those 
of  physical  baselines,  one  would  use  the  inverse  of  the  sum  of  the  variances  of  all  of  the 
antennas  involved. 


5.4  Data  Quality  Check 

Before  continuing  to  the  Doppler  and  chirp  processing,  the  quantity  in  Eq.  (6)  is  tested.  If 
it  is  more  than  10“81,  the  cosine  determination  is  assumed  to  have  failed.  It  is  possible  that 
after  the  first  least  squares  fit  of  £j  on  the  1200 \/2  foot  baseline,  the  inequality  in  Eq.  (6) 
holds,  so  that  a  final  least  squares  fit  fe  is  not  made,  but  that  the  left  hand  side  is  larger 
than  10-81,  so  that  the  direction  cosine  resolution  is  considered  failed. 
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6  Doppler  and  Chirp 


A  512-point  FFT  is  applied  to  an  array  containing  the  cosine  and  sine  of 

2tt(A  <£(n)  —  A  <f>(n  —  1)), 

the  increments  of  the  phase  history  determined  in  Section  4,  weighted  by 

wtf  n  wtf  n_i 
wtf„  +  Wtf„_i  ' 

The  resulting  power  spectrum  is  smoothed  with  3-point  triangular  smoothing: 

Pi  «=  -22 (1 j  +  2.545 Pi  +  Pi+l), 

and  Co  is  chosen  as  the  bin  number  of  the  peak  power  (only  indices  497-512,  and  1-48  are 
checked),  adjusted  somewhat  by  the  strengths  of  the  neighbor  bins.  Then  a  64-point  FFT 
is  applied  to  the  array  consisting  of  the  cosine  and  sine  of 

2tt  (A <j>(n)  -  .bCotl)  , 

weighted  by  wtfn,  and  the  power  spectrum  is  smoothed.  The  bin  number  of  the  peak  is 
chosen  as  D0. 

The  values  Do  and  C0  just  determined  are  used  as  initial  values  of  the  differential  Doppler 
and  chirp  for  a  least  squares  estimate  best  fitting  the  phase  difference  history,  A<f>(n).  The 
calculated  values  are  given  from  the  model  as 

c(n)  =  A  <j>+  Dtn  +  - Ct2n , 

with  A<£  initially  A<^(0),  the  adjusted  phase  difference  profile  value  at  frame  NFTO  from 
Section  4,  and  D  and  C  the  differential  Doppler  and  chirp,  initially  D0  and  CQ.  The  partial 
derivatives  are 

dc(n )  dc(n)  dc(n)  1  2 

dA<f>  ~  1  dD  =  tn  dC  =  2  n' 

These  form  the  nth  row  of  the  matrix  A.  The  normal  equation  is  then 


E-tfn 

(  1  ft) 

tn  *1  \Pn  , 

(AA  <>\ 

|  AD  1  = 

'  £wtfn(A  <f>(n)-c(n))  > 

£wtfnfn(A<f>(n)  -  c(n)) 

U<n  \tl  \tn) 

\  AC) 

^£wtfn^(A<^(n)  -  c(n))y 
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where  the  sums  run  over  n.  The  first  iteration  is  carried  out  without  a  tolerance  requirement 
for  the  residuals.  For  the  second  iteration  only  observations  meeting  the  tolerance 

wtf„(A<£(n)  -  c(n))2  <  T0\ 


are  allowed,  where 

T0i  =  4  max  |  > 

and  S2  is  the  weighted  sum  of  the  squared  residuals  from  the  first  iteration. 


After  the  second  iteration,  if  the  dot  product  of  the  quantities  AlW(0  —  C)  and  the 
applied  correction  is  more  than  20%  of  S2,  the  weighted  sum  of  the  squared  residuals  which 
met  the  tolerance  during  the  second  iteration,  a  third  iteration  is  made,  with  a  new  tolerance: 


Tqi  =  4  max 


/  —  A 

lS„(So-3)’  /’ 


where  So  is  the  number  of  observations  that  met  the  tolerance  in  the  second  iteration. 


The  final  value  of  D  is  multiplied  by  the  sample  rate  and  added  to  the  receiver  bin  center 
frequency  to  yield  the  Doppler.  The  chirp  is  C  times  the  square  of  the  sample  rate. 

The  bin  sizes  for  Doppler  and  chirp  (before  the  triangular  smoothing)  are  0.86  Hz  and 
5.9  Hz/sec,  respectively,  with  ranges  of  ±27.5  Hz  and  [—94.5,284.6]  Hz/sec,  respectively. 
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A 


Least  Squares 


The  term  least  squares  is  used  to  mean  any  process  of  minimizing  a  squared  difference. 
The  best  known  example  is  fitting  a  curve  to  measurements,  but  there  are  many  scientific 
and  engineering  applications  of  the  technique,  which  may  be  viewed  as  a  generalization  of 
the  curve-fitting  problem.  The  theory  is  developed  in  many  places  including  pp.  123-5  of 
[Bate],  and  will  not  be  repeated  here.  We  now  summarize  the  method  (and  standardize  the 
notation): 

We  are  given  a  set  of  parameters  {Pj},  which  are  to  be  improved  or  determined  from 
consideration  of  a  set  of  observations  {0,}.  In  the  usual  case  there  are  more  observations 
than  parameters.  The  model  is  a  functional  relationship  (or  “fit”)  O  =  f(P),  so  form  the 
calculated  values,  Cj  =  fj{Pi),  and  let  A  be  the  matrix  of  partial  derivatives: 

A.  -dJl 

13  dp' 

and  let  the  weight  matrix  W  be  the  diagonal  matrix  whose  elements  are  the  weights  assigned 
to  each  observation.  The  theory  says  that  P  can  be  improved  to  P  +  A P  where  A P  is  the 
solution  of  the  normal  equation 


A'WAAP  =  AlW{0  -C). 


In  terms  of  matrix  elements,  we  have 


dP,  dPj 


weightobs, 


and  the  right  hand  side  is  given  by 

[AlW(0  -  Cj).  =  J^-weightobs(0  -  C). 

obs  1  3 


The  matrix  AKWA  is  known  as  the  inverse  covariance  because  (AW^)'1  gives  the  co- 
variance  matrix  for  the  new  estimates.  In  particular,  its  ith  diagonal  entry  is  the  variance 
(square  of  the  standard  deviation)  of  the  estimate  of  the  zth  parameter. 
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